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Abstract —  An  equation-error  (EE)  method  is  described  for 
estimating  the  three  motion  parameters  of  a  watercraft  (sound 
source)  moving  on  the  sea  surface  at  a  constant  speed  in  a 
constant  direction  as  it  transits  past  a  hydrophone  located  above 
the  sea  floor  in  a  shallow  water  environment,  using  multipath 
delay  measurements  from  the  hydrophone.  The  motion 
parameter  estimates  are  obtained  via  a  linear  least-squares 
(LLS)  minimization  followed  by  an  algebraic  transformation.  A 
weighting  matrix  is  derived  for  the  LLS  minimization  to  reduce 
the  error  variances  of  the  motion  parameter  estimates. 
Computer  simulation  results  show  that  the  error  variances 
approach  the  Cramer-Rao  lower  bounds  for  small  multipath 
delay  measurement  errors.  The  effectiveness  of  the  EE  method  is 
demonstrated  using  real  hydrophone  data.  A  method  is  then 
described  for  estimating  the  linear  trajectory  of  the  transiting 
source  using  the  motion  parameter  estimates  from  the  individual 
hydrophones  of  an  underwater  acoustic  sensor  network. 

I.  Introduction 

Passive  sonar  systems  can  be  used  to  detect  and  localize 
fast  surface  craft  that  generate  intense  continuous  broad-band 
sound  underwater  as  a  result  of  propeller  cavitation  [1],  In  a 
shallow  water  environment,  the  signal  emitted  by  a  surface 
acoustic  source  arrives  at  a  hydrophone  located  above  the  sea 
floor  via  a  direct  path  and  multipaths.  Given  the  sensor  height 
and  water  depth,  the  instantaneous  range  of  the  surface 
acoustic  source  from  the  sensor  can  be  estimated  using  the 
difference  in  time  of  arrival  (or  multipath  delay)  measurement 
between  the  direct  path  and  bottom-reflected  path  signals  [1]- 
[3].  Also,  if  the  source  moves  at  a  constant  speed  in  a  constant 
direction  as  it  transits  past  the  sensor,  then  the  three  source 
motion  parameters  (speed  together  with  time  and  range  at 
which  the  source  is  at  the  closest  point  of  approach  (CPA)  to 
the  sensor)  can  be  estimated  by  observing  the  multipath  delay 
over  a  sufficiently  long  time  period  that  covers  both  inbound 
and  outbound  legs  of  the  source  transit.  A  common  approach 
is  to  fit  a  multipath  delay  model,  which  is  a  nonlinear  function 
of  the  three  motion  parameters,  to  the  sequence  of  multipath 
delay  measurements  in  a  least-squares  sense  [1],  The  three 
motion  parameter  values  that  minimize  the  squared  deviations 
of  the  multipath  delay  observations  from  their  predicted 
values  provide  the  nonlinear  least-squares  (NLS)  estimates  of 
the  speed,  CPA  time  and  CPA  range  of  the  source.  The 
minimization  is  performed  using  numerical  optimization 
methods  which  are  not  only  computationally  intensive  but  also 
require  good  initial  estimates  of  the  three  motion  parameters 
for  fast  convergence  to  the  NLS  solution. 


An  alternative  approach  to  the  problem  is  the  equation- 
error  (EE)  method,  which  has  been  used  previously  to 
estimate  the  depth,  CPA  time  and  CPA  range  of  an 
underwater  acoustic  source,  assuming  the  source  speed  was 
known  [4].  The  main  advantages  of  this  approach  are  that 
closed-form  solutions  exist  (thus  requiring  no  initial  estimates) 
and  it  is  computationally  efficient.  In  this  paper,  an  EE 
method  is  proposed  to  estimate  the  speed,  CPA  time  and  CPA 
range  of  a  surface  craft  (sound  source)  using  multipath  delay 
measurements  from  a  single  hydrophone.  The  method  consists 
of  a  linear  least-squares  (LLS)  minimization  where  a  three- 
dimensional  state  vector,  whose  elements  are  simple  algebraic 
functions  of  the  speed,  CPA  time  and  CPA  range  of  the  source, 
is  estimated.  An  algebraic  transformation  then  converts  the 
state  vector  estimate  into  the  three  motion  parameter  estimates. 
In  order  to  reduce  the  error  variances  of  the  motion  parameter 
estimates,  a  weighting  matrix  is  derived  for  the  LLS 
minimization.  The  error  variances  obtained  from  computer 
simulations  of  the  proposed  method  with,  and  without, 
weighting  are  compared  with  the  Cramer-Rao  lower  bounds 
(CRLBs).  The  proposed  method  is  applied  to  real  hydrophone 
data  and  typical  results  are  presented  to  illustrate  its 
effectiveness.  A  method  is  then  described  for  estimating  the 
linear  trajectory  of  the  transiting  source  using  the  state  vector 
or  motion  parameter  estimates  from  the  individual 
hydrophones  of  an  underwater  acoustic  sensor  network. 

II.  Motion  Parameter  Estimation  -  Single  Sensor 
A.  Problem  Formulation 

Figure  1  shows  the  geometry  of  a  hydrophone  located  at  a 
depth  of  d  below  the  sea  surface  and  a  moving  source  which 
emits  continuous  broadband  sound  underwater  as  it  transits 
past  the  sensor  along  a  linear  trajectory  at  a  constant  speed  v 
on  the  sea  surface.  The  sea  floor  is  assumed  to  be  locally  flat 
around  the  sensor  so  that  the  sea-bottom-reflections  of  the 
sound  emitted  by  the  source  to  the  sensor  during  the  passage 
of  the  source  past  the  sensor  can  be  considered  as  if  they 
originated  from  a  mirror  image  source  moving  at  a  constant 
depth  equal  to  twice  the  local  water  depth  d„ .  The  source  and 
the  sensor  are  located  at  a  height  of  h,  and  above  the  local 
sea  bed  respectively.  The  local  source  height  /;,  equals  dw 
and  the  local  sensor  height  hr  equals  dw-d  .  At  time  rc ,  the 
source  is  at  CPA  to  the  sensor  and  its  horizontal  range  from 
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the  sensor  is  dc .  The  slant  range  of  the  source  from  the  sensor 
at  time  t  is  given  by 

R{t)  =  [V\t-zc)2+R2J2  (1) 

where  Rc  =[(ht  -h,.)2  +d2]l/2  is  the  source’s  slant  range  at 
CPA.  Squaring  and  then  expanding  both  sides  of  (1)  gives 
R2(t)=hT(t)x  (2) 

where  h(t)  =  [t2  ,-2t,\]T  ,  x  =  [v2,rcv2,r2v2 +i?2]r  and  the 

superscript  T  denotes  matrix  or  vector  transpose. 


source 


Fig.l.  Geometry  of  a  sensor  located  above  a  locally  flat  sea  floor  and  passage 
of  a  surface  sound  source  past  this  sensor.  The  image  source  is  located  at  a 
depth  equal  to  twice  the  local  water  depth  or  local  source  height. 


Denote  the  difference  in  the  times  of  arrival  of  the  direct 
path  and  bottom-reflected  path  signals  at  the  sensor  (simply 
referred  to  as  the  multipath  delay)  as  r(t) : 

T(t)  =  [R,.(t)-R(t)]/c  (3) 


where  R,.(t)  =  [R2(t)  +  4h,hr]1/2  is  the  length  of  the  bottom- 
reflected  path  and  c  is  the  speed  of  sound  propagation  in  water. 
It  can  be  shown  using  (3)  that 


m 


4h,hr  -c2z2(t) 
2  cz{t) 


(4) 


Let  i(t)  be  an  estimate  of  the  multipath  delay  z (t)  at  time  t. 
Given  the  sensor  height  hr,  source  height  /?,  and  sound  speed 
in  water  c,  an  estimate  of  R2(t)  at  time  t,  denoted 
as  y(t)  =  R2(t),  can  be  obtained  by  substituting  z(t)  for  r (t)  in 
(4)  and  then  squaring  the  result.  Suppose  there  is  an  error  of 
A  z(t)  in  z(t) ,  which  results  in  an  error  of  e(t)  =  A R2(t)  in  y(l) . 
Using  (2),  the  observation  y(t)  at  time  tk  can  be  written  as 

y(tk)  =  hT(tk)x  +  e(tk)’  1  ^k<K  (5) 

where  K  is  the  number  of  observations.  The  observation  time 
period  q  <t<tk  covers  both  inbound  and  outbound  legs  of  the 
source  transit.  Equation  (5)  can  be  expressed  in  matrix  form  as 
y  =  Hx  +  e  (6) 


where  y  =  [y(tx),...,y(tK)]T  is  the  observation  vector, 
e  =  [e(q ),..., e(tK)]T  is  the  observation  error  (or  equation  error) 
vector,  H  =  [h(q ),...,  h(/g-)]7  is  the  system  matrix,  and 
x  =  [x,,x2,x3)r  is  the  (constant)  state  vector.  The  elements  of 
x  are  simple  algebraic  functions  of  the  three  motion 
parameters  v,  zc  and  Rc ,  and  vice  versa: 

x,=v2,  x2  =  tcv2  ,  x3  =  z2v2  +  R2  (7) 

v  =  x1/2  ,  zc  =  x1~1x2 ,  Rc  =(x3  -xf‘x22)1/2  .  (8) 

The  objective  is  to  estimate  x  for  a  given  y. 

B.  Algorithm 

A  LLS  estimate  of  the  state  vector  x  is  obtained  by  finding 
the  vector  x  that  minimizes  the  following  cost  function: 

J  =  (y  -Hx)r  W(y  -  Hx)  (9) 

where  W  is  a  positive  definite,  symmetric  weighting  matrix. 
The  LLS  estimate  x  is  given  by 

i  =  (HrWH)  HrWy  .  (10) 

Once  x  is  computed,  the  EE  estimates  of  the  three  motion 
parameters,  denoted  as  v,  zc  and  Rc ,  are  readily  obtained  by 
substituting  x  for  x  in  (8). 

If  W  is  an  identity  matrix,  then  x  is  an  unweighted  LLS 
estimate.  Using  a  proper  weighting  matrix  can  reduce  the  error 
variances  in  x  .  If  the  observation  errors  e(tk )  can  be  modelled 
as  zero-mean  random  variables  with  a  covariance  matrix  Re , 
a  proper  weighting  matrix  will  be  the  inverse  of  Re .  For  small 
multipath  delay  measurement  errors  A r(tk) ,  the  observation 
error  e(tk )  =  AR2  (tk )  can  be  approximated  to  the  first  order  as 
e(4 )  =  2  R(tk  ){dR/d  r)  A  r(tk )  (11) 

where  the  derivative  dR/dz  can  be  computed  using  (4)  as 
dR/dz  =  -cRr(tk)/[Rr(tk)~R(tk)\  ■  (12) 

Substituting  (12)  into  (11)  gives 

e(tk)  =  -2 c[R (4 ) - R;1  (4 XT1  A r(4 ) .  (13) 

Assuming  that  all  multipath  delay  measurement  errors  Ar(4 ) , 
1  <k<K  ,  are  independent,  zero-mean,  with  a  common 
variance  of  a2  ,  then  the  observation  errors  e(tk)  are  also 
independent,  (approximately)  zero-mean,  and  the  covariance 
matrix  Re  is  a  diagonal  matrix  whose  kkt h  element,  which 
equals  the  variance  of  e(tk ) ,  is  given  by 

Re(k,k)  =  4c2[R-'(tk)-R;l(tk)]-2a2  fori  <k<K.  (14) 

With  W  =  Rj\  it  can  be  shown  using  (6)  and  (10)  that  the 
mean-square-error  (MSE)  of  x  ,  defined  as  the  ensemble 
average  of  (x-x)(x-x)7  ,  is  equal  to  (F^R^HE1 .  In  practice, 
Re(k,k)  is  computed  using  (14)  with  R(tk)  and  Rr (4) 
replaced  by  their  respective  estimated  values  y1/2(4)  and 
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[y(tk)  +  4hthr]'/2  .  Also,  the  constant  4 c2cr;  in  (14)  can  be 
ignored  as  it  has  no  effect  on  x  . 

C.  Computer  Simulations 

In  the  first  scenario,  the  sensor  height  hr  =  1  m,  source 
height  (equal  to  water  depth)  h,  =  20  m,  speed  of  sound  in 
water  c  =  1520  m/s,  source  speed  v  =  10  m/s,  CPA  time  rc  =  0s, 
and  CPA  slant  range  Rc  =  50  m.  The  multipath  delay  r(t)  was 
computed  using  (1)  and  (3)  every  0.5  s  over  the  time  interval 
of  -10<t<10s.  Multipath  delay  measurements  were  then 
generated  by  adding  independent  zero-mean  Gaussian  noise 
with  a  standard  deviation  of  aT  to  r (?) .  For  a  given  value  of 

(7 T ,  1,000  simulation  runs  were  performed,  first  with  W  =  R ,, 1 
(weighted  EE)  and  then  with  W  equal  to  an  identity  matrix 
(unweighted  EE),  and  the  bias  errors,  standard  deviations 
(STDs)  and  root-mean-square  errors  (RMSEs)  of  v,  z::  and  Rc 
were  computed.  Table  1  shows  the  simulation  results  obtained 
with,  and  without,  weighting  for  (a)  crr  =  10  ps  and  (b) 
<jt  =50  ps.  Also  included  in  Table  1  are  the  results  obtained 
using  the  single-sensor  NLS  method  [1]  and  the  CRLBs.  Both 
the  EE  and  NLS  methods  were  implemented  in  MATLAB®. 
The  NLS  method  used  the  EE  estimates  to  initialize  the 
numerical  (iterative)  minimization.  For  the  number  of 
multipath  delay  measurements  taken  (41  in  each  simulation 
run),  the  computation  time  for  the  NLS  method  was  about  100 
times  longer  than  that  for  the  EE  method. 

Table  1 .  Simulation  results  for  EE  and  NLS  methods  using  single 

SENSOR,  TOGETHER  WITH  CRLBS  ON  STDS,  FOR  TWO  DIFFERENT  MULTIPATH 
DELAY  MEASUREMENT  ERRORS:  (A)  <T r  =  10|jSAND(B)  Gt  =  50  |iS. 


q 

ii 

o 

T: 

Unweighted  EE 

Weighted  EE 

NLS 

CRLB 

V  (m/s) 

bias 

0.030 

-0.070 

0.004 

STD 

0.178 

0.115 

0.114 

0.114 

RMSE 

0.181 

0.135 

0.114 

C  (s) 

bias 

-0.001 

0.000 

0.001 

STD 

0.094 

0.053 

0.053 

0.052 

RMSE 

0.094 

0.053 

0.053 

Rc  (m) 

bias 

-0.048 

-0.050 

0.001 

STD 

0.725 

0.341 

0.340 

0.329 

RMSE 

0.727 

0.344 

0.340 

(b)  C7r  =  50(j.s 

Unweighted  EE 

Weighted  EE 

NLS 

CRLB 

V  (m/s) 

bias 

0.910 

-1.276 

0.017 

STD 

1.245 

0.617 

0.607 

0.569 

RMSE 

1.542 

1.417 

0.607 

A  (s) 

bias 

-0.008 

-0.001 

0.001 

STD 

0.536 

0.342 

0.269 

0.261 

RMSE 

0.536 

0.342 

0.269 

Rc  (m) 

bias 

-1.916 

-1.455 

0.001 

STD 

7.762 

1.653 

1.648 

1.645 

RMSE 

7.994 

2.202 

1.648 

Simulations  have  also  been  performed  for  a  second 
scenario  where  Rc  was  decreased  to  25  m  while  other 
conditions  remained  unchanged.  Similar  trends  to  Table  1 
have  been  observed.  The  computer  simulation  results  show 
that  the  NLS  estimates  are  better  than  the  EE  estimates;  not 
only  their  bias  errors  are  smaller  but  their  STDs  are  also  closer 
to  the  CRLBs.  This  is  not  surprising  as  the  multipath  delay 
measurements  contained  additive  independent  zero-mean 


Gaussian  noise  and  so  the  NLS  estimates  are  also  maximum 
likelihood  estimates.  The  STDs  of  the  weighted  EE  estimates 
approach  the  CRLBs  as  the  value  of  <rr  decreases,  and  their 
accuracy  is  better  than  that  of  the  unweighted  EE  estimates, 
especially  for  Rc  when  the  value  of  <rr  is  large. 

D.  Experimental  Results 

In  a  shallow  water  experiment  (water  depth  ~  20  m),  eight 
hydrophones  were  located  at  a  nominal  height  of  1  m  above 
the  sea  floor.  The  sensor  configuration  was  (almost)  linear 
with  an  intersensor  spacing  of  14  m.  A  small  surface  vessel 
travelled  in  a  straight  line  at  a  nominal  speed  of  9  knots  (4.63 
m/s)  past  the  hydrophone  array  at  a  horizontal  distance  of  less 
than  90  m  from  the  centre  of  the  array.  The  vessel  trajectory 
was  apparently  perpendicular  to  the  array  axis.  The  output  of 
each  hydrophone  was  sampled  at  250  kHz.  The  data  recorded 
from  each  hydrophone  were  processed  in  non-overlapped 
blocks,  each  containing  131,072  samples  (about  0.52  s  of 
data).  Each  data  block  was  subdivided  into  five  75% 
overlapped  blocks,  each  consisting  of  65,536  samples.  The 
power  spectra  of  these  five  subdivided  data  blocks  were 
computed  using  the  Fast  Fourier  Transform  (FFT)  and  then 
averaged  to  produce  a  smoothed  power  spectrum.  An  inverse 
FFT  was  applied  to  the  logarithm  of  the  smoothed  power 
spectrum  to  produce  a  cepstrum.  The  time  lag  at  which  the 
cepstrum  attains  its  maximum  value  provides  an  estimate  of 
the  multipath  delay.  In  this  way,  a  sequence  of  multipath  delay 
estimates  was  obtained  for  each  hydrophone. 


Fig.  2.  Experimental  results  showing  the  sequence  of  multipath  delay 
estimates  from  sensor  6  and  the  prediction  computed  using  the  estimated 
motion  parameter  values. 

The  proposed  (weighted)  EE  method  was  applied  in  turn  to 
the  sequence  of  multipath  delay  estimates  from  each  sensor. 
Figure  2  shows  (as  dots)  the  sequence  of  multipath  delay 
estimates  from  sensor  6  and  (as  a  solid  line)  the  prediction 
computed  using  the  estimated  source  motion  parameter  values 
from  that  sensor.  Figure  3  shows  (as  dots)  the  EE  estimates  of 
the  speed,  CPA  time  and  CPA  slant  range  of  the  surface  vessel 
from  each  of  the  eight  sensors.  Also  included  (as  circles)  in 
Fig.  3  are  the  three  motion  parameter  estimates  obtained  using 
the  single-sensor  NLS  method  for  comparison  purposes.  The 
two  sets  of  estimates  are  in  good  agreement.  The  speed 
estimates  from  both  methods  also  agree  closely  with  the 
nominal  speed.  Note  that  since  the  vessel  trajectory  is 
apparently  perpendicular  (possibly  at  a  small  inclination  angle) 
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to  the  array  axis,  its  CPA  horizontal  ranges  to  any  two 
adjacent  sensors  should  differ  by  a  value  equal  approximately 
to  the  intersensor  spacing,  i.e.,  14  m.  The  CPA  slant  range 
estimates  obtained  by  both  methods  were  converted  into  CPA 
horizontal  range  estimates  and  then  the  difference  in  the 
estimates  for  every  two  adjacent  sensors  was  computed.  For 
the  EE  method,  the  mean  and  standard  deviation  in  these 
differences  are  14.11  and  1.02  m  respectively,  while  for  the 
NLS  method,  they  are  14.10  and  1.01  m  respectively. 


Fig.  3.  Experimental  results  showing  both  EE  and  NLS  estimates  of  (a)  speed, 
(b)  CPA  time  and  (c)  CPA  range  of  the  surface  vessel. 


III.  Motion  Parameter  Estimation  -  Multiple  Sensors 


A.  Source-Sensor  Model 

Figure  4  shows  the  general  geometrical  configuration  for  a 
network  of  N  widely  distributed  underwater  acoustic  sensors 
and  a  moving  source  which  emits  continuous  broadband 
sound  underwater  as  it  transits  over  the  sensor  network  along 
a  linear  trajectory  at  a  constant  velocity  V  on  the  sea  surface. 
The  XT-plane  coincides  with  the  planar  sea  surface.  The 
position  of  sensor  n  is  given  by  ( Xn,Yll,d„ ) ,  where  dn  is  the 
depth  of  sensor  n,  for  1  <  n  <  N  .  It  is  assumed  that  sensor  1  is 
located  directly  below  the  origin  so  that  Xt  =Yt  =0  .  The 
trajectory  of  the  source  on  the  XT-plane  is  described  by 


X(l)  =  dc  cos 0C  +  ((-  tc)V sin  6C 
Y(t)=dc  sin<?c  -(t-rc)Vcos0c 

where  rc  is  the  time  when  the  source  is  at  CPA  to  sensor  1; 
dc  and  0C  are  the  respective  horizontal  range  and  azimuth 
angle  of  the  source  at  CPA  to  sensor  1.  The  source  trajectory 
is  specified  by  the  four  motion  parameters:  V ,  tc,  dc  and  0, . 
Note  that  V  can  be  negative;  the  sign  convention  is  that  V  is 
positive  (negative)  when  sensor  1  is  on  the  right  (left)  hand 
side  of  the  source  as  the  source  moves  along  its  trajectory. 


CPA  at  tc 


Fig.  4.  General  geometrical  configuration  for  a  network  of  N  underwater 
acoustic  sensors  and  transit  of  a  surface  sound  source  over  the  sensor  network. 

The  sea  floor  is  assumed  to  be  locally  flat  around  each 
sensor  -  see  Fig.  1.  At  time  rc  „ ,  the  source  is  at  CPA  to  sensor 

n  and  its  horizontal  and  slant  ranges  from  this  sensor  are  drn 
and  Rc  „  respectively,  for  1  <  n  <  N  ,  with  rcl  =  tc  ,  dc  l  =  dc 
and  Rcn  =  R,  .  The  source  and  sensor  n  are  located  at  a  height 
of  htn  and  hr  n  above  the  local  sea  bed  respectively,  where 
htn  equals  the  local  water  depth  dw  „  and  hl%n  =  „  -  dn . 


B.  Algorithm 

The  EE  method  described  in  the  previous  section  can  be 
applied  to  the  sequence  of  multipath  delay  measurements  from 
each  sensor  to  estimate  the  state  vector  and  hence  the  speed  as 
well  as  the  time  and  slant  range  at  which  the  source  is  at  CPA 
to  that  sensor.  The  sequence  of  multipath  delay  measurements 
from  each  sensor  should  be  taken  during  the  passage  of  the 
source  past  that  sensor.  Thus  the  measurement  times  (i.e., 
instants  at  which  the  multipath  delay  measurements  are  taken) 
and  number  of  measurements  may  vary  from  sensor  to  sensor. 
Only  the  source  speed  \V\  can  be  estimated  with  a  single 
sensor.  Let  v  =|  V  |  and  denote  the  state  vector  for  sensor  n  as 

x„  =K •  Parallel  to  (7)  and  (8),  xln,x2<n  and  xx„ 
are  related  to  v,  rcn  and  Rc  „  via 

xi,„  =  v2  >  *2.„  =  A,„v2  ,  xx„  =  r2„v2  +  i?2„  (16) 

V  =  XTn  •  Tc,n  =  XUnX2,n’  Rc,n  =  (*3,„  ~  •  (17) 

Since  R]n  =  rf2„  +(h,  „  „  )2 ,  it  follows  by  substituting  the 

expression  forT?^^  in  (17)  into  this  expression  that 

dc,n  =[*3,n  ~xlnx2,n  ~(ht,n  ~hr,n)2^2  ■  (18) 
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Denote  the  (weighted)  LLS  estimate  of  the  state  vector  x„ 
from  sensor  n  as  x„  =[xu,,x2n,xXn]T  and  the  corresponding 
EE  estimates  of  the  source  speed  v,  CPA  time  rcn  and  CPA 
horizontal  range  dcn  as  vn  ,  icn  and  dcn  respectively,  for 
1  <  n  <  N  .  These  EE  estimates  are  computed  by  substituting 
x„  for  x„  in  (17)  and  (18).  An  improved  estimate  of  v  is 
obtained  by  averaging  the  source  speed  estimates  from  the 
individual  sensors,  i.e.,  v  =  N~iYN ,  v„  .  The  estimates  of  the 
CPA  time  rc ,  denoted  as  ic ,  and  CPA  horizontal  range  dc , 
denoted  as  dc ,  are  simply  given  by  rcl  and  dc  X  respectively. 

An  EE  approach  is  adopted  to  obtain  an  estimate  of  the 
CPA  azimuth  angle  0C  using  the  state  vector  or  source  motion 
parameter  estimates  from  the  N  sensors.  Consider  the  source 
trajectory  and  the  projection  of  the  sensors  onto  the  XT-plane 
as  shown  in  Fig.  5,  where  sensor  1  is  located  at  the  origin  and 
sensor  n  at  (Xn ,  Y„ ) .  The  unit  directional  vector  v  is  defined 

as  v  =  [sin  9C  -  cos  0C  ] T  so  that  the  source’s  travel  direction  is 
given  by  sgn(F)v  ,  where  sgn(F)  is  the  sign  of  V  .  The  unit 
directional  vector  u  =  [cost?,., sin 0C]T  is  obtained  by  rotating  v 

by  90°  in  the  anti-clockwise  direction.  From  Fig.  5,  it  can  be 
shown  for  2 <n<N  that 

dc,\  -dc,n  =rju  (19) 

V(rc,„-Tc,  i)  =  rjv  (20) 

where  r„  =  [Xn,Y„]r  is  the  position  vector  of  sensor  n. 


Fig.  5.  Source  trajectory  and  projection  of  sensors  onto  the  XT-plane. 


Substituting  V  =  sgn(K)v  into  (20)  gives 

V(.Tc,n~Tc,  i)=rJ«  (21) 

where  the  unit  vector  a  =  sgn(K)v  indicates  the  source’s  travel 
direction.  Using  (17),  the  left  hand  side  of  (21)  can  be 
expressed  as 

v(^,„-rc,i)  =  /(x„)-/(Xl)  (22) 

where  f(x)  =  xxi/2x2.  Replacing  {v,rc  „}  and  {v,  zc} }  by  their 
estimates:  {vn,icn}  and  {vj , f c a }  modifies  (21)  to 

vJc,n-Viic.  i=rZa  +  wn-u  2<n<N  (23) 

where  wn_x  is  the  error  in  observing  v(rc  n  -  zc  ] )  (or  equation 
error).  Equation  (23)  can  be  written  in  matrix  form: 


z  =  Sa  +  w  (24) 

where  the  vectors  z  =  [v2ic2  -Vjfc ,!,•••, %A,iv  ~V\tc,i]T  and 
w  =[wl,...,wN_l]T  ,  and  the  matrix  S  =  [r2,...,rA,]r  .  The  LLS 
estimate  of  a  is  given  by 

d=(SrW„S)-1SrW„z  (25) 

where  W„  is  a  weighting  matrix.  Once  a  =  [a, ,  a2  ] T  is 
computed,  the  unit  vector  u  can  be  estimated  using  the 
relations  between  u  and  v  and  between  v  and  a  as  u  =  [il , ,  ?7  2  ] T  , 
where  ux  =-sgn (V)a2  and  u2  =sgn(K)a1  ,  provided  that  the 
sign  of  V  is  known.  The  estimate  of  the  CPA  azimuth  angle, 
denoted  as  9C ,  is  then  given  by  0C  =  tan  1  (w2 /wj ) . 

If  the  observation  errors  w„_x  can  be  modelled  as  zero- 
mean  random  variables  with  a  covariance  matrix  Rw  ,  a 
proper  weighting  matrix  for  (25)  will  be  the  inverse  of  R,„ . 
The  observation  error  w„  ,  in  (23)  is  caused  by  the  errors  in 
the  state  vector  estimates  x,  and  xn  .  Using  (22),  it  can  be 
approximated  (to  the  first  order)  as 

w„_l=VTf(i„) Ax„  -Vr/(ii)A*,  (26) 

where  V/'(  x  t )  is  the  gradient  of  /(x)  evaluated  at  x/f,  and 
Axt  is  the  error  in  xk,  for  1  <k<N.  By  definition, 
Axk  =xk-xk.  For  small  multipath  delay  measurement  errors, 
a  first-order  approximation  gives  £[Axt.]»0  and  hence 
E[wn_  i  J  =  0  .  Therefore,  R„,  *£[ww!]  .  Assuming  that  Axk 
and  Ax,  are  independent  for  k^l  and  1  <k,l<N ,  it  can  be 
shown  using  (26)  that 

E[wm_xwn_x}  =  WT /{i^ElAxM  ]V/(  i,)  (27) 

E[wlx]  =E[wm_xw^]  +  Vr/(xn  )£[Ax„Ax^  ]V/(i„)  (28) 

for  m*n  and  2  <m,n<N .  In  (27)  and  (28),  .EjAx^Ax^]  is 
the  MSE  of  x^  ,  which  is  equal  to  (Hj.R^H,.  )_1 ,  where  Ht 
and  Rf  t  are  the  system  matrix  and  covariance  matrix  of  the 

observation  error  vector  e  for  sensor  k,  respectively. 

The  sign  of  the  source  velocity  V  is  yet  to  be  determined. 
Substituting  dc  l  ,  dc„  and  u  for  dc  l  ,  dc  ll  and  u  in  (19) 
results  in  an  equation  error  given  by 

=  <4,i  ~dc,n  -sgn(F)rJ [~a2,al]T  .  (29) 

Expressing  (29)  for  2 <n<N  in  matrix  form  gives 

£(y)  =  b-sgn(F)S[-a2,a,]r  (30) 

where  the  vectors  b  =  [dcl-dca,...,dcl-dcN]T  and 
e(V)=[sl(V),...,sN_l(V)]1  .  Define  the  cost  (or  error)  function 
p(V)  =  et (V)e(V)  and  denote  the  estimate  of  f  as  V  .  If 
p(y)  >  p(-v) ,  then  V  =  -v  (V  <  0) ;  otherwise  V  =  v  (V  >  0) . 
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C.  Computer  Simulations 

Five  sensors  were  located  at  a  depth  of  19  m  below  the  sea 
surface  in  a  cross  shaped  configuration.  The  (X,T)-coordinates 
of  sensors  1  to  5  were  (0,0),  (-50,0),  (50,0),  (0,50)  and  (0,-50) 
m  respectively.  A  surface  acoustic  source  travelled  through 
the  network  of  sensors  and  the  values  of  the  source  motion 
parameters  were  K  =  10m/s,  zc  =0s,  dc=  25  m  and  9C  =90°. 
The  sea  floor  was  flat  throughout  the  area  where  the  sensors 
were  located,  and  the  water  depth  was  20  m.  For  each  sensor, 
a  noisy  multipath  delay  measurement  was  generated  every  0.5 
s  over  a  20  s  time  interval  centred  at  the  time  of  CPA  to  that 
sensor,  with  the  measurement  errors  being  independent,  zero- 
mean  Gaussian  distributed  with  a  common  variance  of  <rr . 
For  a  given  value  of  <jt  ,  1,000  simulation  runs  were 
performed  and  the  statistics  including  bias  errors,  STDs  and 
RMSEs  of  V,  ic,dc  and  6C  were  compiled. 

Figure  6  shows  the  series  of  multipath  delay  measurements 
from  each  sensor  in  a  typical  simulation  run  for  o .  =  1 0  gs. 
Table  2  shows  the  compiled  statistics  of  the  source  motion 
parameter  estimates  for  (a)  aT=  10  us  and  (b)  <rr  =  50  gs. 

The  RMSE  of  6C  increases  with  the  value  of  aT  from  0.34° 

to  2°.  Also  included  in  Table  2  are  the  statistics  obtained  using 
the  NLS  method  [5]  and  the  CRLBs  for  comparison  purposes. 
The  NLS  method  used  the  source  motion  parameter  estimates 
from  the  proposed  method  as  the  initial  estimates.  Though  the 
NLS  method  provides  more  accurate  estimates  than  the 
proposed  method,  the  latter  provides  closed-form  solutions 
and  is  computationally  more  efficient  than  the  former. 
Computer  simulations  have  also  been  performed  for  another 
scenario  where  9C=  45°  while  other  conditions  remained 
unchanged.  Similar  results  to  Table  2  were  obtained. 
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Fig.  6.  Sequence  of  multipath  delay  measurements  from  each  sensor  in  a 
typical  simulation  run  for  crr  =  10  ps. 


IV.  Conclusions 

Multipath  propagation  can  be  utilized  for  localizing  fast 
surface  watercraft  in  a  shallow  water  environment.  An  EE 
method  has  been  described  for  estimating  the  speed,  CPA  time 
and  CPA  range  of  a  sound  source  moving  on  the  sea  surface  at 
a  constant  speed  in  a  constant  direction  as  it  transits  past  a 
hydrophone  located  above  the  sea  floor,  using  multipath  delay 


measurements  from  the  hydrophone.  The  performance  of  the 
EE  method  was  studied  by  computer  simulations.  Using  a 
proper  weighting  matrix  reduces  the  error  variances  and 
improves  the  accuracy  of  the  EE  estimates.  Experimental 
results  were  presented  to  demonstrate  the  effectiveness  of  the 
EE  method.  A  method  has  also  been  described  for  estimating 
the  four  motion  parameters  (which  specify  the  linear  trajectory) 
of  the  transiting  source  using  the  EE  estimates  from  the 
individual  hydrophones  of  an  underwater  acoustic  sensor 
network.  Computer  simulations  produced  encouraging  results. 
The  main  advantages  of  the  proposed  methods  over  the  NLS 
methods  are  that  closed-form  expressions  exist  for  the  motion 
parameter  estimates  and  they  are  computationally  efficient. 
The  proposed  methods  can  be  used  to  provide  initial  estimates 
for  the  NLS  methods.  Future  work  will  include  studying  the 
effects  of  the  sensor  configuration  on  the  performance  of  the 
proposed  method  for  a  network  of  sensors  and  verifying  the 
effectiveness  of  the  method  using  real  hydrophone  data. 

Table  2.  Simulation  results  for  EE  and  NLS  methods  using  five 

SENSORS,  TOGETHER  WITH  CRLBS  ON  STDS,  FOR  TWO  DIFFERENT  MULTIPATH 
DELAY  MEASUREMENT  ERRORS:  (A)  C7r  =  10|iSAND(B)  <Tr  =  50  gS. 


(a)  <JT  =  10|ts 

EE 

NLS 

CRLB 

V  (m/s) 

bias 

-0.052 

-0.000 

STD 

0.053 

0.023 

0.023 

RMSE 

0.075 

0.023 

A  09 

bias 

0.001 

-0.000 

STD 

0.025 

0.013 

0.013 

RMSE 

0.025 

0.013 

dc  (m) 

bias 

0.002 

0.001 

STD 

0.207 

0.091 

0.091 

RMSE 

0.207 

0.091 

0C  (rad) 

bias 

0.000 

0.000 

STD 

0.006 

0.002 

0.002 

RMSE 

0.006 

0.002 

(b)  <J T  =  50gs 

EE 

NLS 

CRLB 

V  (m/s) 

bias 

-1.090 

0.009 

STD 

0.292 

0.121 

0.117 

RMSE 

1.128 

0.121 

TC(S) 

bias 

0.000 

-0.003 

STD 

0.142 

0.064 

0.064 

RMSE 

0.142 

0.064 

dc  (m) 

bias 

-0.136 

-0.008 

STD 

1.069 

0.459 

0.456 

RMSE 

1.078 

0.459 

&C  (rad) 

bias 

-0.001 

0.000 

STD 

0.035 

0.012 

0.012 

RMSE 

0.035 

0.012 
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